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If dark matter (DM) is composed by particles which are non-gravitationally coupled to ordinary 
matter, their annihilations or decays in cosmic structures can result in detectable radiation. We 
show that the most powerful technique to detect a particle DM signal outside the Local Group is 
to study the angular cross-correlation of non-gravitational signals with low-redshift gravitational 
probes. This method allows to enhance signal-to-noise from the regions of the Universe where the 
DM-induced emission is preferentially generated. We demonstrate the power of this approach by 
focusing on GeV-TeV DM and on the recent cross-correlation analysis between the 2MASS galaxy 
catalogue and the Fermi-LAT gamma-ray maps. We show that this technique is more sensitive than 
other extragalactic gamma-ray probes, such as the energy spectrum and angular autocorrelation 
of the extragalactic background, and emission from clusters of galaxies. Intriguingly, we find that 
the measured cross-correlation can be well fitted by a DM component, with thermal annihilation 
cross section and mass between 10 and 100 GeV, depending on the small-scale DM properties and 
gamma-ray production mechanism. This solicits further data collection and dedicated analyses. 


I. INTRODUCTION 

The origin of cosmic structures is well understood in 
terms of evolution of matter perturbations arising after 
the inflationary period. Inhomogeneities starting off with 
higher-than-average density grow through gravitational 
instability. DM is a necessary ingredient to the process, 
as it provides the potential-wells where standard mat¬ 
ter is accreted after decoupling and protohalos form. As 
structure formation evolves, DM halos of increasing size 
form in a bottom-up fashion. 

If DM is in form of particles which exhibit non- 
gravitational couplings to ordinary matter, a certain level 
of emitted radiation is expected. Photons can be pro¬ 
duced from interactions of DM with the ambient medium 
(e.g., through scatterings) or from DM annihilation or 
decay by means of direct emission or through the pro¬ 
duction of intermediate particles. The non-gravitational 
signal associated to decay is proportional to the DM den¬ 
sity: it is stronger at low redshift, because the produced 
radiation is diluted by the expansion of the Universe more 
rapidly than its source, i.e. the DM particle density. The 
DM annihilation signal, which is proportional to the den¬ 
sity squared, is also peaked at low redshift since the den¬ 
sity contrast associated to cosmic structures grows non- 
linearly. 
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DM constitutes the backbone of all cosmic structures 
and DM halos represent, collectively, a potential source 
of DM decay or annihilation signals. This means that 
even if the radiation originating from DM annihilations 
or decays in a single halo is too faint to be detected, their 
cumulative signal and its spatial coherence could be. In 
addition, since the DM signal is expected to peak at z < 
0.3, it can be separated by more mundane astrophysical 
processes that typically trace the star formation history 
and peak at higher redshifts. 

To increase the sensitivity to non-gravitational DM 
sources one needs to isolate the annihilation/decay sig¬ 
nal produced at low redshift. An effective way to filter 
out any signal that is not associated to DM-dominated 
structures or that is originated at high redshift is to cross¬ 
correlate the radiation held with bona fide low-redshift 
DM tracers CHS]. In the following, we adopt this ap¬ 
proach in the specific and yet very relevant framework of 
weakly interacting massive particles (WIMP) that may 
either annihilate or decay. We will use the results of 
the cross-correlation analysis between 7 -ray maps from 
Fermi-LAT [?) and the 2MASS catalogue of relatively 
nearby galaxies [5] presented in 0 . 


II. DATA AND MODELS 

The cross angular power spectrum (CAPS) between 
the unresolved 7 -ray sky observed by Fermi-LAT and the 



2 


distribution of 2MASS galaxies can be written as [3]: 

cf 9) = I g Wj (x) Wgix) p 19 (k = e/x, x) , (1) 

J X. 

where x( z ) denotes the radial comoving distance, 
W, (x) represent the window functions described below, 
P-,g(k, z) is the three-dimensional cross power spectrum 
(PS), k is the modulus of the wavenumber, and £ is the 
multipole. Indices 7 and g refer to 7 -ray emitters and 
extragalactic sources in 2MASS, respectively. In Eq. ([l]) 
we used the Limber approximation m , since P ig varies 
(relatively) slowly with k. 

The (differential in energy) window function for 7 -ray 
emission from DM annihilation W^{z) is [3]: 


W?(z) = 


(flDMPc) 2 {< JgV ) 
87T 771dm 2 


(1 + z ) A 2 ( 2 :) 


dN a 
d E-y 
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where S 2 dm is the daQ mean density in units of the 
critical density p c , A 2 (z) is the clumping factor, ttidm 
is the mass of the DM particles, and (o a v) denotes the 
velocity-averaged annihilation rate. dN a /dEj indicates 
the number of photons produced per annihilation and 
determines the 7 -ray energy spectrum. The exponential 
damping quantifies the absorption due to extra-galactic 
background light [ 12 ] . 

The window function for DM decay is [3J: 


w dl z \ = ^DMpcPd d N d c _ t[<|EtW1 
7 477 777dm d Ey 


(3) 


where Td = I /77 is the DM decay rate. 

The window function of 2MASS galaxies is W g (z) = 
H(z)/cdN g /dz and their redshift distribution dN g /dz is 

m 


dN g p z m 

^ (z)= r(f)r exp 



(4) 


with 777 = 1.90, p = 1.75 and Zq = 0.07. 

We employed the 2MASS catalogue instead of other 
compilations because the galaxy distribution in Eq. [4] is 
peaked at very low redshift as for the DM emission of 
Eqs. [2] and [3] This enhances the cross-correlation signal. 
The picture for astrophysical components would be dif¬ 
ferent and other catalogues might be more informative 
(see Section S2 in the Supplemental Material can- 

The PS P ig in Eq. (JT|) is computed within the halo- 
model framework, as the sum of a one-halo plus a two- 
halo terms. For more details, see [§]. Both the PS and 
the clumping factor A 2 (z) in Eq. j2l) depend on a num¬ 
ber of DM properties: the halo mass function, that we 
take from Ref. [22] , the halo density profile, for which we 


1 A 6-parameter flat ACDM cosmological model is assumed with 
the value of the parameters taken from Ref. eh. 


assume a Navarro-Frenk-White model I2.'{l . the minimum 
halo mass, that we set equal to 10~ 6 Mq, and the halo 
mass-concentration relation c{M,z ), that we adopt from 
Ref. [23j- The theoretical uncertainty of these quanti¬ 
ties is rather small for halos larger than 10 10 Mq, be¬ 
cause they can be constrained by observations and simu¬ 
lations. Since the DM decay signal is mainly contributed 
by large structures, the theoretical predictions are rel¬ 
atively robust. This is not the case for the annihilation 
signal which is preferentially produced in small halos and 
in substructures within large halos. Consequently, theo¬ 
retical uncertainties on the annihilation signal are larger. 
For the subhalo contribution we consider two scenarios 
(low and high) to bracket theoretical uncertainty. The 
LOW case follows the model of Ref. [25] (see their Eq. (2), 
with a subhalo mass function dn/dM su b oc M“^). The 
high scenario is taken from Ref. 26j, with the halo mass- 
concentration relation extrapolated down to low masses 
as a power law. Further uncertainties in the concentra¬ 
tion parameter and in the value of the minimum halo 
mass can introduce an extra factor of ~ 2 of uncertainty, 
that we quantify only in the Supplemental Material m , 
for the sake of definiteness. 

In our CAPS model (Eq. IT]) , we add a constant term 
Cih ( one-halo correction term) to correct for possible un¬ 
accounted correlations at very small-scales, within the 
Fermi-LAT Point Spread Function (PSF). The value of 
C\h will be determined by fitting the data, and we an¬ 
ticipate that we find a C\h value compatible with zero. 
Thus, the inclusion of this term does not change signif¬ 
icantly the results. For a discussion on this term, see 
Ref. [5]. 

The measured CAPS C^ 19 ^ is a convolution of the true 
CAPS and the effective beam window function Wf that 
accounts for the PSF of the instrument and the pixeliza- 
tion of the 7 -ray map. Both quantities depend on en¬ 
ergy. We use the Wf derived in Ref. [9] and model the 
observed spectrum as cj, 19 ' 1 = Wf C^ 19 ' 1 . 

In the following, we shall consider the angular cross 
correlation function (CCF) rather than the spectrum. To 
model the CCF, we Legendre-transform the CAPS: 

CCF^ 9 \e) = Y ^-^C7 9 P([ cos(0)] , (5) 

t 

where 9 is the angular separation and Pe are the Legendre 
polynomials. 

To compare model and observed CCFs, we estimate 
the x 2 difference defined as: 

X 2 = EE K - m oM)) [Clei ~ m l^) ’ 

n—1 6i 0j 

, (6) 

where m and d indicate model and data, n identifies 
each one of the three overlapping energy ranges consid¬ 
ered (E > 0.5,1, and 10 GeV) and the indices 9 t and 
9j run over 10 angular bins logarithmically spaced be¬ 
tween 9 = 0.1° and 100°. C g g . is the covariance matrix 
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FIG. 1. Cross-correlation above 500 MeV for the best fit¬ 
ting annihilating and decaying DM scenarios, compared to 
the measured CCF. The curves are for DM particles of 100 
GeV (200 GeV) annihilating (decaying) into bb. We show 
the two annihilation models HIGH and LOW with annihi¬ 
lation rates (a a v) = 2 x 10~ 26 cm 3 s _1 (blue-dashed) and 
2.4 x 1CF 25 cm 3 s _1 (blue-solid), respectively, and a decay 
model with lifetime r = 1.6 x 10 27 s (red-dotted). The green 
curve shows the CCF of the 1-halo correction term C\h- We 
show the sum of this component and the DM CCF (in the 
LOW scenario) with the black curve. The inset shows that 
these DM models provide a subdominant contribution to the 
observed IGRB spectrum |27j . 


that quantifies the errors of the data and their covariance 
among the angular bins. Data and covariance matrix are 
taken from Ref. [9]. The parameter vector for annihi¬ 
lating DM is A = [todm, (vav), C\h], whereas for the 
decaying DM is A = [?n D M, r d , Ci h }- 


III. RESULTS 

In Fig. [l] we show a comparison between the measured 
CCF in one of the considered energy bins (E > 500 MeV) 
and the best fitting annihilating and decaying DM mod¬ 
els obtained from the analysis discussed below. Error- 
bars are given by the diagonal elements of the covari¬ 
ance matrix. DM models fit the measured CCF remark¬ 
ably well (for best fitting model, Xbf = 16.7 with 26 
d.o.f.). It is also noteworthy that the level of annihi¬ 
lation/decay rate provides a minor contribution to the 
Isotropic Gamma Ray Background (IGRB) measured by 
the Fermi-LAT E3, as shown in the inset of the figure. 
This implies that the cross-correlation technique can de¬ 


tect a DM signals too faint to show up in the total inten¬ 
sity measurement (for a review of the IGRB properties, 
see Ref. pS]). 

In Fig. [2j we show the ler and 2er C.L. contours (ob¬ 
tained marginalizing over C\h) for DM mass and anni¬ 
hilation/decay rate for various final states. Note that, 
although we use only three energy bins, they are suffi¬ 
cient to constrain the DM mass which induces a small 
but characteristic signature in the energy spectrum. In 
the LOW scenario the ler region lies just above the thermal 
annihilation rate ( a a v) = 3 x 10 ~ 26 cm 3 s _1 . In the high 
case, the DM signal increases by a factor of ~ 10 and con¬ 
sequently regions shift down by one order of magnitude. 
Therefore, given the current uncertainty in modeling DM 
structures we conclude that the thermal cross section is 
well within the allowed regions for uidm 200 GeV. 

We stress that the confidence contours in Fig. [2] are 
drawn under the assumption of no contribution from as- 
trophysical sources. While their purpose is mainly il¬ 
lustrative, they may not be unrealistic since astrophysi- 
cal sources, that are indeed required to account for the 
IGRB thanks to their medium-to-large redshift emission, 
can indeed provide a negligible contribution to the cross¬ 
correlation signal between Fermi-LAT and 2MASS galax¬ 
ies that, as we point out, has a rather local origin (see dis¬ 
cussion in Section S2 of the Supplemental Material [HD- 
On the other hand, given the current uncertainty on the 
astrophysical components of the IGRB, an astrophysical 
model that can explain the measured cross correlation 
signal with no additional contribution from DM can be 
found |9] . Future data and analyses will help distinguish¬ 
ing between these two options. 

This cross-correlation measurement can alternatively 
be used to derive 95% C.L. upper bounds on the an¬ 
nihilation/decay rate. These bounds are conservative 
and robust, since we assume here that DM is the only 
source of the 7 -ray signal, without introducing additional 
assumptions on astrophysical components which would 
make make the constraints stronger but also more model 
dependent. The 95% C.L. upper bounds on the WIMP 
annihilation (decay) rate as a function of WIMP mass 
are shown in the left-hand (right-hand) panel of Fig. [3] 
For bb and t + t~ final states, the thermal annihilation 
rate is excluded for masses below 10 (100) GeV in the 
LOW (high) scenario. In the case of the bounds 

degrade by about one order of magnitude. 

In Fig. [4] we compare the sensitivity of our cross cor¬ 
relation method with that of other extragalactic 7 -ray 
probes. We focus on these probes since they are similarly 
affected by uncertainties in modeling DM halo and sub 
halos properties. This allows to compare various tech¬ 
niques in a homogeneous and robust way, something that 
cannot be done with local DM tracers (Galactic regions, 
dwarf galaxies) or early Universe probes, which have dif¬ 
ferent systematic uncertainties (see however the discus¬ 
sion in Section SI of the Supplemental Material [14]). For 
illustrative purposes, we selected the LOW substructure 
scheme and bb final states case. We verified that differ- 
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FIG. 2. Left: la and 2a allowed regions for the DM annihilation rate versus its mass, for different 7 -ray production channels 
and assuming a LOW substructure scheme. Crosses indicate the best fitting models. In the HIGH scenario, regions remain similar 
but shifted downward by a factor of ~ 12, see stars indicating the best fitting models. Right: The same as in the left-hand 
panel but for decaying DM, showing the DM particle lifetime as a function of its mass. 




FIG. 3. Left: 95% C.L. upper limits on the DM annihilation rate as a function of its mass. Both HIGH and LOW clustering 
schemes are shown for WIMPs annihilating into bb (with the impact of different energy bins reported for the latter case). Other 
final states of annihilation (p + fj,~ , t + t ~ , W + W~) are shown in the LOW scenario only, for clarity. Right: 95% C.L. lower 
limits on the DM lifetime as a function of its mass, for different final states of decay. 


ent choices provide little differences and the results are 
robust to both the DM clustering model and the anni¬ 
hilation/decay channel. We consider again the simplest 
case (where most conservative bounds can be derived), in 
which the astrophysical contribution is set to zero in all 
observables and only DM is contributing as 7 -ray source. 

The bound corresponding to the IGRB energy spec¬ 
trum has been derived using the IGRB estimated by the 
Fermi-LAT Collaboration [27] and adding up in quadra¬ 
ture statistical and systematic errors given in their Table 
3. For the autocorrelation bound, we considered the an¬ 


gular spectrum estimated in four energy bins in Ref. [29] 
as provided in their Table II (DATA:CLANED) and av¬ 
eraged in the multipole range 155 < l < 504. For both 
probes, the model prediction has been computed using 
the same DM modeling as in our analysis. Our bounds 
are compatible with the ones presented in Refs. [5UH55] 
(under the same set of assumptions). Cluster bounds are 
instead taken directly from the literature. In particu¬ 
lar, for annihilating DM, we consider the analysis of 34 
clusters using expected sensitivity for the 5 years Fermi- 
LAT data in Ref. [36] which uses the same LOW model 
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FIG. 4. Ratio £ = {a a v) i /(a a v) x (in the annihilating DM 
case) and £ = tx /t; (in the decaying DM case) between the 
95% C.L. bounds derived with method discussed in this work 
(cross-correlation of 7 -rays with 2MASS catalogue, labeled 
with X) and from other extragalactic 7 -ray probes (i stands 
for total IGRB intensity, angular autocorrelation or clusters). 
The plot refers to the bb final state and (for the annihilating 
DM case) the low substructure scheme. 

adopted here. For decaying DM, we consider the anal¬ 
ysis of 8 clusters in 3 years of Fermi-LAT data taking 
performed by [37j. 

Fig. i shows that the cross-correlation technique 
stands out as the most sensitive one, improving the con¬ 
straints by a factor between a few to a hundred over the 
other techniques. Note that the ratio decreases at high 
energy because our analysis focuses at low energies (up 
to E > 10 GeV). Since the IGRB is measured up to 820 
GeV there is room for further improvements. 

IV. CONCLUSIONS 

We compared the predicted angular cross-correlation 
between the 7 -ray emission induced by DM annihilation 
or decay and the distribution of 2MASS galaxies with the 
measured CCF between these objects and the Fermi-LAT 
7 -ray maps. 

The contribution of astrophysical sources to the IGRB 


is assumed to be subdominant at low redshift and not in¬ 
cluded in the model prediction, so to derive conservative 
bounds on DM. We found that in the LOW [25] and high 
[25] scenarios the “thermal” annihilation cross-section is 
excluded at 95% C.L. up to DM masses of 10 and 100 
GeV, respectively, for final state of annihilation into bb 
and t + t ~ . 

We demonstrated that the cross-correlation technique 
is significantly more sensitive to a DM signal than all 
other extragalactic 7 -ray probes used so far. This was 
done by comparing the bounds of our cross-correlation 
analysis with the most recent results from IGRB, angular 
autocorrelation and clusters, finding an improvement of a 
factor ranging from a few up to 100 for both annihilating 
and decaying DM. 

We showed that a WIMP DM contribution can fully 
explain the observed cross correlation. A canonical 
WIMP with a mass in the 10-100 GeV range, annihi¬ 
lation rate around the thermal value, and realistic model 
for DM halo and sub-halo properties reproduce both 
size and shape of the measured angular cross-correlation. 
This intriguing possibility deserves further investigation 
within a more comprehensive framework that include 
contributions from astrophysical sources and additional 
data. 

Future investigation employing the Pass 8 release 
from the Fermi-LAT and forthcoming surveys at 
low/intermediate redshift [35] will therefore provide re¬ 
markable insights to the particle DM quest. 
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Supplemental Material 


SI. COMPARISON WITH OTHER 7-RAY 
BOUNDS 


When discussing bounds on WIMP DM, one typically 
refers to the constraints on its microscopic properties, 
i.e., on the mass and annihilation (decay) rate, for a given 
set of branching ratios of annihilation (decay) into ordi¬ 
nary particles. Especially in the case of indirect signals, 
such bounds strongly depend on the DM distribution on 
different astrophysical scales: this often suffers of sizable 
uncertainties, which in turn reflect into uncertainties on 
the DM bounds. For this reason, when comparing the 
results of this work with other bounds from the litera¬ 
ture, we restricted our discussion in the main text to ex- 
tragalactic probes: in fact, the bounds arising from the 
7 -ray average intensity, from the 7 -ray angular autocor¬ 
relation and from galaxy clusters, have an origin which is 
(at least partially) common to the cross-correlations sig¬ 
nal we discuss here, and therefore their magnitudes are 
expected to scale similarly when assuming different DM 
models. For this reason we have included them explicitly 
in the main text. 

In this Section, we attempt instead to compare our re¬ 
sults with bounds that originates from 7 -rays produced 
within or nearby our own Galaxy, focussing on annihilat¬ 
ing DM for definiteness. Since this comparison is more 
model dependent and much less straightforward, we de¬ 
fine, for each probe, a conservative (CONS) and opti¬ 
mistic (OPT) scenario of DM spatial distribution (paral¬ 
leling the LOW and high scenarios of the main analysis). 
We stress that this discussion only aims at providing a 
generic (although informative) comparison, whilst a com¬ 
prehensive assessment of the various uncertainties in each 
single channel is clearly beyond the scope of the present 
analysis. 

As mentioned, in the case of the angular cross¬ 
correlation between 7 -rays from Fermi-LAT and 2MASS 
galaxies, the OPT and CONS scenarios are the HIGH and 
LOW cases introduced in the main text. The associated 
95% C.L. bounds are reported with solid lines in Fig. SI 
(taken from Fig. 3) for a bb final state. 

Bounds from 7 -ray emission in the halo of the Milky 
Way (MW) in Fig. SI are taken from Ref. [15]. To quan¬ 
tify uncertainties, we consider two different values of the 
local DM density, namely po = 0.7 GeV cm -3 (OPT) and 
Po = 0.2 GeV cm " 2 3 (CONS).^ 

Dwarf spheroidal (dSph) satellites of the MW are prob¬ 
ably the objects that currently set most stringent 7 -ray 
limits on DM. The red dashed line in Fig. [Si] is taken 


2 The reported curves correspond to the 3cr bounds for an NFW 
profile, as shown in Fig. 4 of Ref. CD- Considering the 95% 
C.L. limits or a different choice for the DM profile would not 

significantly alter the curves shown in Fig. IS 11 


from the analysis of Ref. mi. which includes ultra¬ 
faint dSphs. Although the addition of the recently dis¬ 
covered ultra-faint dSphs can improve the DM bounds, 
the low number of detected stars in these systems makes 
the derivation of their dispersion velocity profile (and in 
turn of the DM content) quite uncertain. Moreover, their 
proximity to the MW suggests strong tidal interactions 
which are normally not accounted for when solving the 
Jeans equation at the equilibrium to derive the DM pro¬ 
file. A more conservative assessment of the uncertainties 
associated to the DM profile of ultra-faint dSph galaxies 
makes them less constraining than classical dSphs (see, 
e.g., Ref. [TTj [TH])- We thus include only the classical 
ones in the CONS scenario (and, conservatively, assume 
a Burkert profile which reduces the bounds by a factor 
of ~ 1.3, see Fig. 7 of Ref. iIB]). In order to re-derive 
the Fermi-LAT bounds for this scenario, we rescale the 
bound of Ref. |T 6 ] by the ratio of the largest J-factors be¬ 
tween the ultra-faint and classical dSphs. Although this 
rescaling is simplistic, nevertheless it grabs the essence 
of the uncertainty. We also note here that other (typi¬ 
cally ignored) uncertainties can affect the J-factor of the 
same classical dSphs by a factor of a few (e.g., system- 
atics from triaxiality, see Ref. El) : this implies that the 
blue dashed bound in Fig. [Si] could be further pushed 
upward. 


Finally, we consider the limits coming from the Galac¬ 
tic Center (GC). They strongly depend on the assump¬ 
tion concerning the DM density in the innermost region 
of the MW, something which is not currently obtainable 
through observations. This fact makes the bounds associ¬ 
ated to the GC rather uncertain. The dot-dashed curves 
shown in Fig. SI are taken from Ref. m which makes 
use of four years of Fermi-LAT data. The OPT scenario 
refers to a compressed NFW, while the CONS case to 
a canonical NFW. Note that they span nearly three or¬ 
ders of magnitude. For comparison, we show also the 
preferred region for the so-called GC excess EU j. 

As for the DM bounds derived in the main text, and 
reproduced in Fig. SI we adopted the case where only 
DM is considered to contribute to the 7 -rays signal. In 
this way, the derived bounds are conservative, and no 
further assumption on the various possible astrophysical 
backgrounds (which would make the comparison more 


model dependent) is required. As clear from Fig. SI 


the bounds derived in this work are competitive with 
the most stringent Galactic bounds. This also implies 
that a DM explanation of the cross-correlation signal (see 
Figs. 1 and 2 in the main text), or more generally DM 
scenarios lying in the region between the high and LOW 
constraints from cross-correlation, are not currently ruled 
out by other means. 

We finally comment that we restricted the discussion to 
7 -ray bounds only. We have not introduced other probes 
(such as charged cosmic-rays, radiative emissions, neu¬ 
trinos, CMB) since they require further and completely 
different sets of assumptions (e.g., on the astrophysi¬ 
cal environment and galactic transport) and one can al- 









FIG. SI. Bounds on the DM annihilation rate as a function 
of its mass for various 7 -ray probes, for a bb final state. The 
limits derived in this work from extragalactic DM signals are 
compared with the constraints from dSph galaxies [T5], the 
MW halo [15] and the GC [15]. For all cases we define an 
optimistic (OPT) and conservative (CONS) scenario. The 
region compatible with the so-called GC excess | 20 | is shown 
for comparison. 

ways select a reasonable scenario such that the associated 
bounds become weaker than the ones shown in Fig. [ST] 
for the cross-correlation analysis. 

In conclusion, the limits derived in the main text from 
the analysis of the cross-correlation signal between 7 -rays 
extragalactic emission and large scale structure tracers 
are competitive and independent from the other bounds 
that can be derived from alternative means, and, as 
shown in the main text, represent the most constraining 
ones when compared to bounds arising from extragalactic 
probes. 

S 2 . DISCUSSION ON 7-RAY ASTROPHYSICAL 
BACKGROUNDS AT LOW 2 

The extragalactic 7 -ray signal from annihilating or de¬ 
caying DM peaks at very low redshift (in Fig. S2, we show 
the contribution to the IGRB as a function of z for the 
same annihilating DM LOW case of Fig. 1; see, e.g., the 
discussion on DM window functions in Ref. [3] for more 
details). This implies that the size of the signal of its 
cross-correlation with a low redshift catalogue is tightly 
connected to the total DM contribution to the IGRB. 
This is not the case for more mundane astrophysical 7 - 
rays sources, as already mentioned in the Introduction of 


the main text. Indeed, the fraction of IGRB emitted by 
astrophysical sources at z < 0.1 is below few percent of 
their total IGRB contribution, as can be clearly seen in 
Fig. S2 (the models shown in Fig. S2 are the same em¬ 
ployed in Ref. [5T] and fit the whole Fermi-LAT IGRB 
measurement 1221). This is a general property that holds 
for all the main extragalactic 7 -ray source populations. 

Therefore the tight link between cross-correlation sig¬ 
nal and contribution to the IGRB that holds true for 
DM is not valid, at low redshifts, for the astrophysical 
sources. A large contribution to the IGRB might be ac¬ 
companied by a low cross-correlation signal at low red¬ 
shift, and a scenario where astrophysical sources make up 
the whole IGRB (shown in the inset of Fig. 1) with, on 
the other hand, a DM contribution dominating the cross¬ 
correlation with 2MASS (more in general, with a cata¬ 
logue which peaks at very low redshifts, see e.g. Fig. S2) 
is a viable situation. Two questions are at hand here: 
i) whether such a scenario is reasonable from the point 
of view of properties of astrophysical 7 -rays sources; ii) 
whether this options is compatible with data. Let us 
start discussing the first point. 

The vast majority of extragalactic sources present in 
the Fermi-LAT catalogue are blazars. This implies they 
are the 7 -ray emitters whose average properties are best 
known. As shown in Fig. S2, the unresolved emission 
from BL Lacertae (BL Lac) and flat spectrum radio 
quasars (FSRQ) is strongly suppressed at low redshift. 
Therefore blazars cannot explain the measured cross¬ 
correlation of Fermi-LAT data with 2MASS. 

Two other extragalactic populations are thought to 
provide a large contribution to the IGRB: misaligned 
AGN (mAGN) and star forming galaxies (SFG). Unfor¬ 
tunately, since a very low number of these objects have 
been detected so far by the Fermi-LAT telescope, the 
knowledge of their 7 -ray luminosity function suffers of 
large uncertainties, even in the low redshift regime. We 
cannot obtain firm predictions for the expected overall 
size of the associated cross-correlation signal. 

Concerning the angular behaviour, let us note that, at 
scales < 1°, the DM interpretation shown in Fig. 1 stems 
from the non-linear clustering. At the redshift of 2MASS 
(z ~ 0.07) such angular scales correspond to physical 
scales of ~Mpc. Therefore, a component providing a 
good fit to the data needs to have a power spectrum with 
a sizable non-linear term at Mpc scales, which disfavours 
an interpretation in terms of galactic objects. 

In summary, there are currently significant uncertain¬ 
ties on the low-z 7 -ray luminosity function of the astro- 
physical components of the IGRB which prevent a ro¬ 
bust determination of their contribution to the measured 
cross-correlation of Fermi-LAT data with 2MASS. 

In the present analysis, we have been investigating the 
assumption of a subdominant astrophysical contribution 
combined with a particle DM explanation of the cross¬ 
correlation signal. Confronting this scenario with data in 
the main section of the paper, the question ii) have been 
addressed. We have demonstrated that it is, in principle, 
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FIG. S2. Contribution to the IGRB for E 7 > 1 GeV from var¬ 
ious astrophysical 7 -ray emitters as a function of redshift. The 
curves refer to the unresolved component (we set the threshold 
for source detection to F aens = 5 x 1CD 10 photons cm -2 s -1 ). 
The black solid curve shows instead the redshift distribution 
of objects in the 2MASS catalogue. 


a viable possibility. 


S 3 . UNCERTAINTIES FROM THE CHOICE OF 

Mmin AND c(M) 


For the sake of definiteness, in the main text, we fix 
the minimum halo mass to be M m i n = 10 -6 Mq, and 
take the halo mass-concentration relation c(M, z) from 
the reference case in Ref. [55] at masses above 10 1 o Mq. 
The latter is then extrapolated to lower halo masses with 
two different prescriptions, high and LOW, as explained 
above. 

In this Section, we explore the impact of a different 
choice for M m i n and of a different concentration model. 
In Fig. S3, we show the clumping factor A 2 introduced 
in Eq. (2). It is a “normalization” of the DM signal 
which encodes the dependence on DM clustering prop¬ 
erties. The gray band shows the A 2 region obtained by 
varying the c(M ) relation within the band provided in 
Ref. [55]. Dotted lines are instead for M m i n = 10 _ 12 Mq 
( blue) and M min = 1O _ 2 M 0 (red) (the fact they are close 
the borders of the gray band is just coincidental). In this 
case, we focus on the LOW clustering scenario, but similar 
scalings occur for the HIGH case as well. 

Both effects introduce an extra factor of ~ 2 of uncer¬ 
tainty in the predicted signal for annihilating DM. Such 


FIG. S3. Clumping factor A 2 as a function of redshift 
for different choices of minimal halo mass and mass- 

concentration relation c(M), as described in the labels. 


uncertainty band would translate into a corresponding 
uncertainty in ( <j a v) of Fig. SI (solid curves). 

In the case of decaying DM, the normalization of the 
signal is instead not affected by the assumptions on M m i n 
and c(M). 

Above we just discussed the impact on the overall 
normalization of the signal. Varying the c(M) relation, 
the relative contributions of different halo masses would 
change and in turn the angular shape of the correlation 
can be affected. This is however a minor effect (for rea¬ 
sonable c(M)) and is not relevant given the current ac¬ 
curacy of data. 







